function S_crtbp = I2C_u(S_i,m1,m2,DU,thetao,t_s)
%% INPUTS
% S_crtbp = State in CRTBP frame at time t
% t = time in DU;
% m1 = Mass of planet 1 (Central planet)
% m2 = mass of the planet 2 (Moon)
% DU = distance unit
% Wcm = angular speed of the moon
% thetao = intial angular offset from positive x axis at t = 0

%% ASSUMPTIONS
% 1. Moons are in co-planar circular orbit.
% 2. they are all happy.

%% CONSTANTS
G           = 6.674e-11 * ((1/1000)^3);
%[] Gravitational Constant

GM_1 = m1*G;
%[] Gravitational Parameter of first mass

GM_2 = m2*G;
%[] Gravitataional parameter of second mass

LU = DU;
%[] Length unit

TU = sqrt(LU^3/GM_1);
%[] Time unit

VU = LU/TU;
%[] Velocity Unit

Wcm = sqrt((GM_1+GM_2)/LU^3)*[0;0;1];
%[] Angular speed of planet 2 (Moon)

thetadot = sqrt(GM_1*DU)/DU^2;
%[] Theta dot of the planet of interest.
%% Determine position of the moon in Inertial frame.

thetaf = thetao+thetadot*t_s;
%[] Angular offset from thetao for time t

R2_interest = DU*[cos(thetaf);sin(thetaf);0];
%[] Position of the Moon at time t_s

V2_interest = sqrt(GM_1/DU)*[-sin(thetaf);cos(thetaf);0];
%[] Velocity of the moon at time t_s

S2_interest = [R2_interest;V2_interest];
%% Determine center of mass. 
Scm1 = (GM_2*S2_interest)/(GM_2+GM_1);
%[] state of barycenter with respect to primary body at time t_s

Rcm1 = Scm1(1:3);
%[] Position component of Scm1

Ssatcm = S_i-Scm1;
Rsatcm = Ssatcm(1:3);
Vsatcm = Ssatcm(4:6);

%% Determine Teci2crtbp and Tcrtbp2eci

Hm          = cross(R2_interest,V2_interest);
%[1/s]Current specific angular momentum of the Moon WRT the CM in ECI coordinates.

Rhat        = R2_interest / norm(R2_interest);
%[]Current CRTBP radial direction in ECI coordinates.

Nhat        = Hm / norm(Hm);
%[]Current CRTBP normal direction in ECI coordinates.

That        = cross(Nhat,Rhat);
%[]Current CRTBP tangential direction in ECI coordinates.

Tcrtbp2inert  = [Rhat, That, Nhat];
%[]Matrix that transforms vectors from CRTBP to ECI coordinates.

Tinert2crtbp  = transpose(Tcrtbp2inert);
%[] Matrix that transforms vectors from ECI to CRTBP coordinates.


%% Determine S_crtbp

R_crtbp = Tinert2crtbp * Rsatcm;
%[] Determine positon of interest in CRTBP frame.

V_crtbp = Tinert2crtbp * Vsatcm - cross(Wcm,Tinert2crtbp*Rcm1) - cross(Wcm,R_crtbp);
%[] Determine positon of interest in CRTBP frame.

S_crtbp = [R_crtbp/LU;V_crtbp/VU];
%[] State of the S_eci eexpressed at CRTBP frame at JD of interest.


end

